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Abstract 

Slender structures immersed in a cross flow can experi- 
ence vibrations induced by vortex shedding (VIV), which 
cause fatigue damage and other problems. VIV models 
in engineering use today tend to operate in the frequency 
domain. A time domain model would allow to capture 
the chaotic nature of VIV and to model interactions with 
other loads and non-linearities. Such a model was devel- 
oped in the present work: for each cross section, recent 
velocity history is compressed using Laguerre polynomials. 
The compressed information is used to enter an interpola- 
tion function to predict the instantaneous force, allowing 
to step the dynamic analysis. An offshore riser was mod- 
eled in this way: Some analyses provided an unusually fine 
level of realism, while in other analyses, the riser fell into 
an unphysical pattern of vibration. It is concluded that 
the concept is promissing, yet that more work is needed 
to understand orbit stability and related issues, in order 
to further progress towards an engineering tool. 



1 Introduction 

Vortex induced vibration (VIV) is a vibration of an elas- 
tic structure that occurs when a fluid flowing around the 
structure sheds vortices at near-regular intervals, locked 
with the structure's own vibration. VIV is a major con- 
cern in the offshore oil industry in particular, where marine 
currents can cause slender structures like pipelines, risers, 
umbilicals and cables to vibrate, inducing fatigue damage. 
VIV is a hard problem because on one hand full hydrody- 
namic computations of vortex sheddings from structures 
are as yet impractical, and on the oher hand, it is chal- 
lenging to simplify a strongly non-linear dynamic system. 
Semi empirical VIV models provide the state of the art of 
VIV engineering. They work by predicting added mass and 
excitation coefficients on the basis of reduced frequency 
and amplitude of vibration, and seek one or several oscil- 
lation modes that satisfy equilibrium. 

Compared to such semi empirical VIV models, in the long 
term an efficient time-domain VIV model would open new 
possibilities : 

1. Study of VIV on non-linear structures, for example 
studying the damping effect of seafloor interaction in 
a steel riser, or using a hysteretic cross section model 
for VIV on flexible pipes. 



2. Accounting for VIV caused by unsteady water flows, 
in particular by waves or vessel motions. 

3. Accounting for the increase in drag at wave frequency 
due to VIV. 

4. Accounting for the superposition of wave-frequency 
and VIV-frequency stresses in fatigue analysis. 

5. Accounting for the asymmetry of oscillation patterns 
in the vicinity of, for example, a seafloor. 



The objective of the work reported here is to demonstrate 
the viability of a local, deterministic, time-domain force 
model for VIV on slender bodies with cylindric cross sec- 
tions. This force model is used at each Gauss point of 
the dynamic finite element (FE) model of a slender struc- 
ture subject to external steady or unsteady water currents, 
during a time domain analysis (e.g. Newmark-/3 time in- 
tegration with Newton-Raphson iteration). Hence the FE 
model resembles that commonly used in a slender struc- 
ture analysis, with degrees of freedom for the structure, 
and none for the surrounding fluid. In other words, the 
proposed model takes the place usually held in software 
by the Morison model for wave induced loads. 

The litterature describes a few time-domain models of 
VIV that, like the present model, do not explicitely model 
the wake flow. In [8], at any step and point along a 
cable, the recent velocity history is approximated by a 
harmonic function, which is then used to enter charts 
that predict excitation and added mass coefficients as a 
function of reduced amplitude and frequency. A model 
which could be described in the same way, but differs in 
several details was developped by |6J. More recently, the 
hydrodynamic force has been described as the response on 
a non-linear single degree of freedom van der Pol oscilator 
(U [15] [24] [25]. The models enumerated here only deal 
with cross-flow vibration. 



The present model differs from the above ones in that it 
treats in-line and cross flow vibrations jointly, does not 
use a harmonic simplification of the motion or forces, and 
does not reduce the wake response to a single degree of 
freedom system. 



2 Model outline 



2.1 Postulate 

The present work hinges on the following postulate. The 
force exerted by the surrounding fluid on a section of the 
slender structure, is completely determined by the recent 
histories at that section of the velocities of the structure 
and the undisturbed fluid. Several points in this sentence 
are worthy of discussion. 

The "force" includes the components usually distributed 
into added mass, excitation forces, drag, lift etc... . 

That the force "at section of the slender structure" 

is determined by the history "at that section" implies a 
"strip theory" in which it is excluded that motions of the 
structure at a point A cause disturbances in the fluid that 
affects the force at point B away from A. In other words, 
it is assumed that there is no significant transmission of 
information in the axial direction within the water (as op- 
posed to within the slender structure). This would be 
proved wrong if it turned out that unstable phenomena, 
like boundary layer shedding, although transmitting lit- 
tle energy along the structure, transmits information that 
steers how local hydrodynamic energy is channeled at a 
given point along the structure. 

That the force should be "completely determined" im- 
plies that the behavior of the structure is deterministic. 
This does not contradict the observation of hysteretic 
response of short cylinders mounted on elastic support. 
Uniqueness of forces given a position does not imply 
uniqueness of static equilibrium. Neither does "completely 
determined" contradict the observation of irregular and 
unpredictable responses to VIV: non-linear dynamic sys- 
tems can have a chaotic behavior. Still, complete deter- 
minism is provably wrong, since a short vertical cylinder 
dragged at uniform speed through water will experience 
oscillating lift forces. At any given moment, there is noth- 
ing in the history of (constant) velocity that allows to pre- 
dict whether the lift is left or right. So the present work 
is based on the bet that ignoring such "bifurcations" still 
leaves us with a useful model. 

"History" here relates to causality. The force on the 
structure does not depend of future motion of the struc- 
ture. By contrast, frequency-domain models do not make 
an explicit distinction here, typicaly requiring a steady 
state vibration. This can also be contrasted to Morison's 
equation which predicts forces on a cylinder from instan- 
taneous relative velocities and accelerations. 

"Recent" can be defined as anything between the present 
time and a few times t w , where the value of t w still is an 
object of debate. t w is likely to be case dependent. Cur- 
rent will transport (convect) away vortices so that they 
quickly loose significance, so that t w should be of the or- 
der of D/U where D is the cross section diameter and U 



the current velocity. In contrast, if the cylinder is oscil- 
lating in still water, it will be traveling in its own wake, 
and t w should be related to the rate of diffusion and/or 
viscous dissipation of vortices, which is likely to result in 
much higher values of t w . Tests on periodic forced mo- 
tion of short cylinders sometimes show a slow drift of the 
forces (over as many as ten periods). In contrast, force 
decay tests for a cylinder stopped after oscilliations at zero 
mean velocity, point towards a fraction of a period. In the 
present work, the idea is to chose a "universal" value of t w 
for the system, after adequate scaling (cf. Section [372] ) . 

The "velocities" are what count. Accelerations would not 
do because for example, zero acceleration can correspond 
to different speeds and hence different forces. On the 
other hand, the force on a cylinder will not be affected by 
a uniform translation of its whole trajectory, so a history 
of positions contains irrelevant information. 

In the remained of this text, the word "trajectory" will be 
given a very specific meaning. The trajectory is defined 
as the recent history of the velocity vector of the cylinder 
relative to the undisturbed surrounding fluid. 



2.2 Restrictions 

In the present phase of research, the following restrictions 
are introduced, in order to achieve some simplification of 
the task. The outer cross section of the slender structure 
is assumed perfectly circular and smooth. The surround- 
ing fluid is assumed to be infinite, excluding the presence 
of sea floor, free surface or neighbouring risers. Only fluid 
flows perpendicular to the cylinder at any point are con- 
sidered. 



2.3 Input and output 

As stated earlier, the VIV model being developed herel 
replaces the Morison model for wave induced loads. The 
VIV model is called at each step and iteration, and at 
each Gauss point or node of each element. 

The model is to receive as input: 

1. the diameter of the cylinder. 

2. the instantaneous velocity of the cross section rela- 
tive to the undisturbed fluid. 

3. the instantaneous velocity and acceleration of the lo- 
cal undisturbed fluid, in a Galilean reference system. 

The model uses velocity information stored from previous 
steps. On this basis, the model produces as output: 

1. the vector of hydrodynamic forces per unit length, 
acting on the cylinder. 



2. the matrix containing the derivative of the above 
with respect to instantaneous values of the cylinder's 
behaviour. 

Gauss integration is then used to compute a consistent 
load vector and derivative matrix for each element. Note 
that these element matrices are likely to vary significantly 
over each VIV oscillation "period" - in contrast to added 
mass or damping matrices, deemed to be constant over 
a long time in semi-empirical VIV models. The connec- 
tion of the force model to the finite element analysis is 
discussed in Section |7] 



2.4 Algorithmic steps 

Only the local VIV model is described here, not the whole 
FE analysis. 

1. The relative velocity of the cylinder relative to water 
(thereafter: "velocity") is computed. 

2. The velocity is scaled (Reynolds scaling) to that the 
cylinder diameter is the unit of distance (Section 



The trajectory (again: the recent histories of both 
x and y components of velocity) is compressed to a 
small number of "Laguerre coefficients". This com- 
pression is such that it provides accurate information 
over the recent past and increasingly coarse informa- 
tion for more distant past (Section 

The Laguerre coefficients are used to enter an in- 
terpolation function (a feed-forward neural network 
with some specifically tailored properties) which re- 
turns x and y components of hydrodynamic force 
(Section [5J. The fitting of the interpolation function 
is discussed in Section IQI 



5. The force is scaled back to the relevant diameter 



(Section 3.2 ) 



6. The Froude-Krylov forces, which depend on the ac- 
celeration of the undisturbed flow, are added (Sec- 



tion 3.1) 



The identification of non-linear systems using a bank of 
orthogonal filters (including Laguerre filter) to generate 
multiple signals from a single one, and then using the mul- 
tiple signals to enter a non-linear, memory-less function, 
was introduced by Norbert Wiener [26J. In the present 
work, a base of Laguerre polynomials is used, in con- 
trast to Laguerre functions introduced by Wiener. While 
Wiener apparently did not use neural networks as non- 
linear functions (but for example Hermite polynomials), 
neural networks in Wiener models have been studied for 
some time [2J. In the present work, Laguerre filtering is 
presented without making use of the vocabulary of cyber- 
netics. In particular, the z-transform is not introduced 
here. 



3 Ancillary transformations 



3.1 Froude-Krylov forces 

This section gives the justification for point [6] of Sec- 
tion |2.4| If the undisturbed fluid in which the cylinder 
is plunged is accelerating (because of surface waves, for 
example), then it is natural to introduce two reference sys- 
tems: © is a Galilean reference system, for example fixed 
relative to the sea floor, 21 is an accelerated reference sys- 
tem, locally following the undisturbed flow. Transforming 
the equations of equilibrium from (5 (in which we carry 
out FEM analysis) to 21 (for which we have experimental 
data, in water that is no accelerated) requires the addition 
of inertia forces. 

The inertial forces create a uniform pressure gradient that 
was not present in the laboratory test. The effect of a 
pressure gradient on a submerged body is variously re- 
ferred to as "Archimedes forces" when the pressure gradi- 
ent results from the acceleration of gravity, or as "Froude- 
Krylov forces" when the pressure gradient is due to fluid 
acceleration in surface waves. As familiar, the integral of 
the pressure over the wet surface is transformed into a 
volume integral |5j. 

It is assumed that this pressure gradient does not affect 
the turbulent flow, so that the pressure gradient can sim- 
ply be added to the pressures resulting from turbulence. 
This seems reasonable enough for incompressible flows, 
and indeed when it comes to Archimedes forces, the sub- 
merged weight of a cylinder is routinely subtracted to lab- 
oratory measurements and the relevant correction added 
again in FEM analysis - even though the Archimedes 
forces in the laboratory do not necessarily scale with those 
in the analysis. Further there is no experimental indication 
that a horizontal and vertical cylinder, all other conditions 
being equal, experience different forces. 

To conclude, the hydrodynamic force acting on the cylin- 
der at a given instant is the sum of two terms: 



1. A force that is a function of only the cylinder diame- 
ter and the recent history of the velocity of the cylin- 
der relative to the undisturbed, steady water flow. 

2. Froude-Krylov forces. 



All computations in Sections [4] and [5] deal only with the 
first of the above two terms. 



3.2 Scaling 



This section details how points[2]and [5]of Section Z4 are 
implemented. In order to reduce the amount of experi- 
mental data necessary to create the interpolation function 



used in point [4] one must take advantage of scale sim- 
ilarities. To that effect, all data used to either train or 
query the database is scaled. Correspondingly, all forces 
returned by the database are scaled back. 

VIV forces are assumed to be uniquely defined by fluid 
density p, kinematic viscosity v, cylinder diameter D and 
the motion. Hence, in order to create a database that is 
be entered with scaled velocities, we wish all experimental 
data to be scaled to fixed reference values p a , v Q and 
D Q . By expressing the units of these quantities, one gets 
three equations on X m , X s and Xk g , which are the scaling 
factors for the basic units of distance, time and mass. 
Solving the system yields 
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Once the scaling of basic units is known, the scaling of 
any derived quantities e.g. velocities, accelerations and 
forces per unit length can be expressed: 
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Note that since scaling is applied consistently to all 
derived quantities, all non-dimensional numbers based 
on combinations of distance, time and mass (including 
Reynolds and Froude numbers) is conserved. However, 
any dimensional quantity with units different from those 
of p, v and D is scaled to values that depend of p, v and 
D. In particular, Equation |5]shows that all accelerations, 
including the acceleration of gravity g are scaled with a 
factor proportional to D 3 jv 2 . So while the scaling used 
here may conserve Froude's number, it does not allow to 
build a database of forces related to surface wave effects, 
because the database does not refer to a constant value 
9o- 

The choice of p Q , v Q and D Q is arbitrary, and in this work, 
all are set to the value 1. D Q = 1 [m] implies that scaled 
displacements can be considered to have "1 diameter" as 
unit. D — 1 [m] and v Q = 1 [m 2 /s] together imply that 
scaled velocities are expressed as Reynolds numbers since 
the scaled velocity is calculated as Dv/v where wis the 
velocity. 

The Reynolds number is usualy computed using some ve- 
locity characteristic of the system under study. In VIV 
science, the undisturbed velocity of the current is used. 
By contrast, in this work, instantaneous local values of 
the relative velocity vector is multiplied by ~. The scaled 
velocities thus obtained are a generalisation of the tradi- 
tional use of Reynolds number: Considering an immobile 




Figure 1: Weighted Laguerre polynomials (blue) are 
summed (black) to approximate a trajectory (red). Ver- 
tical shifts were added for readability. The 



cylinder in a current, the norm of its scaled relative veloc- 
ity vector is equal to the traditional Reynolds number. To 
prevent confusion of the present usage of Reynolds num- 
ber with the more particular classical one, yet emphasize 
the relation between both, the expression "ilr-Reynolds" 
(for "instant, local, relative Reynolds") will be used in this 
document. 



4 Characterization of trajectory 



4.1 Foreword 



This section details how point [3] in Section 2.4 is to be 
implemented. The objective is, for any given point in 
time, to distill a "summary" of the recent history of the 
velocity of the cylinder relative to the surrounding fluid 
(trajectory). Note that the history of each component of 
the velocity vector is treated separately in this section and 
that the procedure is applied to the scaled trajectory. 

The trajectory is approximated as a linear combination of 
some adequate family of functions, and the coefficients in 
this linear combination are the summary (Figure [I] The 
family of functions that is used here is the series of La- 



guerre polynomials (Sectior 4.2 ). It is shown in Section 
4.3 that if the "Laguerre coefficients" of the linear com- 
bination are obtained by integrating the product of the 
trajectory by adequate "Laguerre dual" functions, then 
the difference between the approximating linear combi- 
nation and the real trajectory is small in the recent past 
and larger in the further past. This justifies the choice 
of Laguerre polynomials: they allow to summarise the 
trajectory in a way that represents recent velocities very 
precisely, and older velocities in a coarser manner. It is 
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were the interpretation of t w has been discussed in Section 




2.1 Functions will now be noted as vectors (in a Hilbert 



space), marked with overlined symbols. An indexed fam- 
ily of functions will be noted as a matrix (symbols with 
double overline) and so will a linear operator (a distri- 
butions of two variables). We introduce the symmetric 
positive definite operator 
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Figure 2: Laguerre polynomials (top), Laguerre duals 
(middle) and weight function (bottom). Taking the con- 
volution of a signal by the Laguerre duals one obtains La- 
guerre coefficients. If one takes the linear combinations 
of the Laguerre polynomials weighted by the coefficients, 
one gets an approximation of the original signal, with a 
quality that decreses towards the past in a way related to 
the weight function. 



assumed that this corresponds to the information needed 
to obtain a good estimate of the hydrodynamic force. 

Computing the integral of the product of Laguerre duals 
and trajectory takes time. Luckily, one can show (Section 
4.5) that the Laguerre coefficients are the solution of a 
differential equation driven by the instant value of the ve- 
locity. To obtain results that are independent of step size, 
this differential equation must be carefully discretized in 



l/l = >/w (12) 

Further we introduce the base 

I(t,i) = d(-t/t w ) t el", ie {1,...,"} (13) 
Equation [8] can be rewritten in matrix notation as 



I = L oWoL 



(14) 



where / is the n x n identity matrix. It is useful to intro- 
duce the weighted norm or iu-norm 
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time (Section 4.6 1 when summarizing experimental data. 



Note that since W(t) is of dimension [l/s], |/| is of 
the same dimension as/. So went taking / as a scaled 
velocity, 1/1 is an ilr-Reynolds number. 



4.3 Analysis and synthesis 



4.2 Definitions 



For a history v(t) of either the x or y component of the 
velocity, we seek the vector of "Laguerre coefficients" f for 
The Laguerre polynomial (Figure [2] top) of degree i - 1 the above base that minimize the weighted discretization 
can be defined by its Rodrigues formula [T| error 

S -») (7) 
Laguerre polynomials verify the orthonormality property 

C t (x) Cj (x) e~ x dx = S i: j (8) 
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We seek to describe the recent trajectory with a precision 
that is good for the immediate past, and decreasing for 
the further past. To this end, we introduce a weight func- 
tion which emphasizes "recent past" (Figure [2] bottom) 



To this effect we require that the derivative be zero 

dJ 

which implies 
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with 



D = W o L (22) 

where the functions D used for analysis consists of the 
Laguerre duals (Figure [2] middle. Not to be confused 
with the Laguerre functions introduced in Equation 27 \ 
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the functions in D and L do not span the same space. 
Hence the appellation "dual base" is abusive. 



4.4 Convergence 



Laguerre functions, which can be defined as 

F(M) = M-t/tw) 
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or, in matrix notation 
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have been extensively studied. Series of Laguerre func- 
tions are known to converge almost everywhere (under 
some conditions of continuity) [18]. In matrix notation 
this result can be stated as 



lim 
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This can be used to obtain a result on the convergence of 
series of Laguerre polynomials. We introduce the change 
of variables 
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We hence have convergence in terms of the quality of 
approximation that we are seeking, with emphasis on the 
recent past. Further, on any finite (or "compact") interval, 
convergence in the w-norm is equivalent to convergence 
almost everywhere. So under some conditions of continu- 
ity on g, the series of Laguerre polynomials obtained using 
D as analysis functions converges almost everywhere to- 
wards g in any finite interval. 

Figure [5] illustrates how Laguerre coefficients indeed pro- 
vide a "summary" of the trajectory 
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Figure 3: Example of Laguerre approximation for two 
components of a velocity history (arbitrary scaling). The 
red dot marks the present time. The red curve is the 
original cyclic signal and the black curves are Laguerre 
approximations for two different instants 

4.5 Differential equation for Laguerre co- 
efficients 

In the finite element analysis, we need to update the La- 
guerre coefficients at each iteration of each time step, for 
every Gauss point of every node of the system. The ex- 



plicit calculation of Equation 21 for every update is hence 
a CPU-time critical operation, taking in the order of nxN 
floating point operations (flops), where n is the number 
of Laguerre polynomial used, and N the number of time 
steps that the dual functions take to decay to a negligible 
value. Further, for each Gauss point, 2N velocity values 
need to be stored, a severe memory requirement. 

In the present Section and the next it is shown how the 



computation of Equation 21 can be carried out by a re 



cursive operation requiring no other storage than that of 
the Laguerre coefficients and the last velocity values, and 
taking in the order of n x n flops, which is advantageous 
because n <C N. In this Section it is shown that r verifies 
a differential equation driven by the history tJof the veloc- 
ity component. In Section 4.6 this differential equation 



is solved time-step by time-step in a recursive update. 



Equation 21 can be rewritten without matrix notation, 
and differentiated 
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A property of Laguerre polynomials is 
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where C i (8) is a generalized Laguerre polynomial 
Hence we can write 
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which is of the form 
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Equation 42 shows that at any time t, the rate of the 



Laguerre coefficients is fully defined by the Laguerre co- 
efficients and the velocity signal. 



4.6 Recursive filter 



The discrete integration of Equation 42 must be done 
carefully, for two reasons. First it is important to obtain 
Laguerre coefficients that are independent of the sampling 
rate used (as long as the sampling rate is "adequate"). 
This is because the experimental data on which the VIV 
model is based may come from experiments which, af- 
ter scaling, may have different sampling rates. Further, 
the numerical analysis in which the VIV model is used 
may use yet another time step. The choice of time step 
or sampling rate must not affect the way a trajectory is 
characterized by Laguerre coefficient. 

The second reason for care in discrete integration is that 
we wish to be able to create synthesized signals L ■ r of 
good quality. Synthesized signal are neither used in the 
numerical process of creating a force interpolation func- 
tion (Section [j>]) or in the FEM use of the VIV model. 



However visualization is essential to the process of re- 
search, both for fault diagnosis and quality control, and 
to communicate an understanding of the method. 

This discrete integration is only used in the analysis of ex- 
perimental data, to provide an input to the training of the 



Votatron" (Section 5.5 1. In dynamic analysis, the integra- 



tion of Equation 42 is done by means of the Newmark-/? 
method, as detailed in Section [7] 

Assume that velocity is sampled at regular intervals 

Vj =v(t + j dt) (45) 

We seek the values of the Laguerre coefficients at the 
same intervals 



r(t +jdt) 
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The vector Tj (the list of the coefficients for all Laguerre 
polynomial, take at step j) must not be confused with 
scalar Tj (the coefficient for the Laguerre polynomial of 
degree i). We choose to such that to + j dt = 0, and we 
approximate v by a function that is linear over the interval 



[0,dt]. Equation 42 becomes 
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This new differential equation can be solved exactly: We 
seek a solution of the form 



: (t) = exp (jit) ■ a + bt + 1 
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over the interval. Here cxp (jit) stands for a matrix ex 
ponential. Replacing this expression into Equation 47 
noting that 
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and identifying the constant and linear terms and enforc- 
ing the initial value leads to 
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Replacing these expressions in Equation 50 at t = dt, a 



tedious but straightforward computation yields the recur- 
sive filter 
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5 Force interpolation 



5.1 Foreword 



This section details the implementation of point[4] in Sec- 
This section presents an interpolation func- 



2.4 



tion 

tion which, given the Laguerre coefficients, predicts the 
present value of the force vector. Polynomials were 
considered initially, but is soon became clear that feed- 
forward "neural networks" provide a better class of func- 
tions to work with. The reason for that is that the number 
of polynomial coefficients of degree d for a polynomial of 
n variables is n d , and high values of d must be expected 
to be necessary. By contrast, in a neural network, non- 
linearity is introduced by "sigmoid" or "threshold" func- 
tions, and the coefficients are used to specify in which 
direction non-linearity applies. Further, polynomials are 
infamous for their propensity to oscillate. 

The "rotatron" presented here is based on the "percep- 
tron" [22 21J, a well studied architecture of neural net- 
work which provides a flexible tool for the interpolation 



of scalar-valued functions of a vector (Section 5.2 1. The 
rotatron takes advantage of certain symmetry properties 



of the physics at hand (Section 5.3) 



In Section [7] the rotatron is used to predict scaled forces 
based on the Laguerre coefficients for scaled trajectories. 



5.2 Perceptron 



A perceptron [22] E] is a simple feed-forward neural net- 
work, consisting of 3 layers. The input layer has 2n neu- 
rons where n is number of Laguerre coefficients for each 
velocity component and the factor 2 comes from the need 
to analyze in-line and cross-flow speed histories together. 
The values of the input layer neurons are set to the La- 
guerre coefficients for both velocity components. The 
second layer has nhid neurons, whose values are an affine 
function of the values of the first layer, passed through a 
sigmoid function like 





Figure 4: Perceptron: "arbitrary" functions (top) can be 
represented as a sum of sigmoid steps (bottom 5) of dif- 
ferent orientation and steepness Njki, shift Vj, and height 



a(x) 



(62) 



15000 




Figure 5: For circular cross section it is assumed that if 
two trajectories can deduced from each other by rotation 
or mirroring, then the corresponding forces are deduced 
from each other by the same operation 



the origin, then the resulting forces are within the same 
line. In particular zero velocities must imply zero forces. 

The symmetries imply that, once experimental data for 
a trajectory has been obtained, there is no need to ac- 
quire data for rotated or mirrored trajectories. However, 
if one was training a perceptron to interpolate the data, 
the training set would need to include trajectories and 
their rotates and mirrors, with the correspondingly ro- 
tated and mirrored forces. This would increase memory 
and CPU usage during training, but also during use of 
the trained perceptron, because the perceptron will need 
a larger number of hidden layer to interpolate the training 
data. 

Another approach is hence used in the present work: the 
classic perceptron is replaced by a "rotatron" (Section 
5.4 1 . It is designed so that, whatever the values of the 



weight coefficient, a rotation or mirroring of the input 
trajectory results in the same rotation or mirroring of the 
output force vector. 



Finally, the third layer gives the output of the perceptron, 
and its values are an affine function of the values of the 
second layer. This can be summarized as 



fi = Mij ■ a (N jk i ■ tm + Vj) + Ui 



(63) 



M^, Njki, Ui and Vj are the "weights" or interpolation 
coefficients, that must be adjusted to fit the perceptron to 
interpolate some given data. Tki are Laguerre coefficients 
and fi are predicted force components, i is the index of 
force direction (x vs. y), j the index of neuron in the 
hidden layer, k the index of velocity direction and I index 
of Laguerre coefficient. 

Each output of the perceptron can be seen as a function, 
which is a sum of sigmoid steps (Figure Q. 



5.3 Symmetries 

The relation between trajectories (in the sense of his- 
tory of the velocity of the cylinder relative to the water) 
and forces can reasonably be assumed to exhibit several 
symmetries (Figure [5]): 

Rotational symmetry If a trajectory can be de- 
duced from the other by a rotation around the origin, 
then the resulting forces are also deduced from each 
other by the same rotation. 

Mirror symmetry If a trajectory can be deduced 
from the other by a mirroring around a line crossing 
the origin, then the resulting forces are also deduced 
from each other by the same mirroring. 

Rotational symmetry and mirror symmetry together, im- 
ply directionality: If a trajectory is within a line crossing 



5.4 Rotatron 



A modified interpolation function (which will be refered to 
as "rotatron" in this text), which enforces the symmetries 
discussed in Section I5.3I is is defined as 



fi = VfcCr 



i/ith 



Oik 

Vjk 



ik 



Dik 



\V[j]k \ (1 + \V[j]k\ ak 



Vlk 



-1 - e 
M M Tji 



Vlk 



(64) 

(65) 
(66) 

(67) 

(68) 
(69) 



In the above, index i and j refer to direction, index I to 
the Laguerre polynomial and k to the hidden layer. Vk, 
Uk and Mki are tunable parameters. Tji are Laguerre co- 
efficients, given as input to the "rotatron". See Appendix 
|A]for conventions on index notations and in particular for 
the syntax |yuife|. Note that Equations ! 
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64 



and 



69 



op- 



erate linearly, identicaly and independently on the terms 
related to the x and y directions, while Equation 66 in- 



volves a unit vector multiplied by a non linear function 
of its norm. Figure [6] illustrates the flow of information, 
from right to left, from two vectors containing the histo- 
ries of the velocity components, to Laguerre coefficient, 
that are then processed in the rotatron. 

The non-linear function appearing in Equation [66] is a 
sigmoid, whose abruptness is parametrized by Uk (Figure 
[7]). The sigmoid is shown in Figure [7J for various values 
of the parameter Uk- 



Figure 6: Laguerre analysis and rotatron transform velocity histories into a hydrodynamic force. The matrix D is the 
discrete form of the Laguerre "duals", which appear in Equation 21 




X 



Figure 7: Log-logistic sigmoid functions 
5.5 Training 

"Training" of a neural network refers to finding weight co- 
efficients Vk, Mfej and, Uk such that for any training point 
number m, consisting of Laguerre coefficients Tji m and 
two force components f im , the outputs fi m computed by 
the neural network are close to fi m . 



5.5.1 Regularization 

A common problem when training neural networks is 
"overspecialization" [23J. In this situation, the neural net- 
work predicts the training outputs with high accuracy but 
behaves wildly between the training points. In contrast, 
what is implicitly sought is a smooth response of the net- 
work to the input, even if this means an imperfect fit to 
the training data. 



regularization [23J: the value of the weight parameters Vk, 
Mki and, Uk are chosen by minimizing the cost function 

J \V[k] ,U[k], M[ fc ] , /[„ n ] , t y Jro ] ) = 

\ (fim ~ fi (r mm )) 2 + P \ {Ul + Vi + M 2 kl ) (70) 

where p is the regularization coefficient, an arbitrary in- 
put to the training algorithm. High values of p favor 
smoothness of the response of the neural network against 
precision in reproducing the training set. 

5.5.2 Conjugate gradient optimization 

J is a function of a large number of weight coefficients, 
and hence it is not practical to compute the Hessian of 
J, because the Hessian is a full matrix. It also proves 
to be very costly to even compute an approximation to 
it as done in the Levenberg-Marquardt algorithm [7] 114], 
On the other hand, the Nelder-Mead "downhill simplex" 
algorithm [19J, which uses only the values of J, proved 
very slow in this case. Hence a search method is chosen, 
that determines the search direction from the gradient of 
J |16| . This is a conjugate gradient method, in which 
the step length is found by deriving the gradient in the 
direction of the search. In this method, the positive def- 
initeness of the (implicit) Hessian is forced by adding a 
scaled identity matrix to it, a technique known as "trust 
region". 

The conjugate gradient method proved far more efficient 
than the Levenberg-Marquardt and Nelder-Mead methods 
for the present task. 



6 Metric 

6.1 Euclidean metric and distance 



Many strategies are described in the literature to address In order to describe the available data, it is useful to de- 
this problem. One of them, which is adopted here, is fine a distance between trajectories. This will allow to 
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determine to what extend the set of available data "fills" 
the set of all possible trajectories, or to detect zones of 
transition from one hydrodynamic behavior to the other. 
Finally, this will help detecting contradictions in the avail- 
able data, arising from a variety of sources, including hid- 
den experimental variables, measurement uncertainties or 
inadequate modeling in inverse methods and not least, 
the natural variability of VIV forces. 

The x and y components of a trajectory are described by 
a pair of functions: 



(71) 



We can define a scalar product between trajectories, that 
captures any recent differences: 




4000 6000 8000 10000 12000 14000 16000 18000 
x-velocity [Re] 



f og= f T x oW°g x + fy oW°g y (72) 



(f x (t)g x (t) + f v (t)g y (t))dt (73) 



By replacing f x , f , g x and g by their expression in 
terms of Laguerre polynomials and their respective La- 
guerre coefficients Tf x , Tf y , r gx and r gv , one finds that 



with 



1 T - 
f 09 



—T — 
T fx ' T gx 

—T — 



Tfx 

T fy 



—T — 

' T fy ' T av 



' gx 



(74) 
(75) 

(76) 



The distance is defined from the scalar product in the 
usual manner: 



\f-g\ = 



\/(f-9) T o(f-g) (77) 
\rf-T a \ (78) 



In other words, neighboring vectors of Laguerre coeffi- 
cients describe trajectories that are similar in the recent 
past. This is illustrated by taking random samples of La- 
guerre coefficients around a given value obtained from 
data analysis and plotting the synthesized trajectories 
(Figure [8]). 



Figure 8: A set of neighboring trajectories according 
to Equation [78] Typical distance between trajectories: 

10 3 [ilrRe]. 

where £H is the set of all rotations of the trajectories 
around the origin and & the set of all mirroring of tra- 
jectories around a line passing by the origin. Note that 
no norm or scalar product associated to the distance d is 
presented here (The vector-space of trajectories, divided 
by the group of rotations and mirrorings, is not a vector 
space). 



Because f x is related to tj x by the same linear relation 
that relates / to Tf y , linear combinations of f x and 
f y (including rotation and mirroring) are related to the 
same linear combinations on Tf x and ft v . By expressing 
the distance / — R(g)\ as a function of the angle a of 
the rotation R, and then differentiating with respect to 
a, it can be shown that the value of a that minimizes 

\7-R(g)\ is 



a = arctan 

( T fx ' T gy 



r fy 



T fy 



Tfx 



(80) 



where arctan (y, x) G ]— 7r.7r] is the angle of a vector 
[x, y] with the cc-axis. Similarly, it can be shown that the 
mirroring that minimizes / — S (g) \ is the composition of 
a rotation of angle 



6.2 Rotatron-distance 



The above does not account for rotational and mirror 
symmetries. We seek a distance for which the distance of 
a trajectory to its transforms by rotation or mirroring is 
zero. Another distance is hence introduced: 



d(f,g) 



min min f 



R{g)\ 



min f 

see 1 



S(g) (79) 



P = arctan 

{ T fx ' Tgy + Tfy • Tg X , 



T fy 



Tfx 



(81) 



by a swap of the sign of the x-coordinates. Equations [80] 



and rH] allow to compute 79 



Figure [9] shows a trajectory and the trajectories within 
a small database that have the smallest distance to it, 
measured using d. 



6.3 Fractal dimension 

Considering a relatively uniform cloud of points, the num- 
ber to of points in a sphere is proportional to the radius r 
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x 10 4 Typical distance 11496 [Re] 



Fractal dimension of p: 4.76 
Fractal dimension of pf: 4.67 




1.5 2 
V x [Re] 



Figure 9: A trajectory and its neighbors in terms of rota- 
tron distance. Smooth curve: Laguerre approximation of 
trajectory, stippled arrow: true force, smooth arrow: pre- 
dicted force. Black is for the trajectory used to enter the 
model to find the force. Green and red are used for the 
three closest points in the database, respectively before 
and after rotation or mirroring. 

of the sphere to the power of p, where p is the dimension 
of the space in which the cloud is defined. For example, 
using 2 x 10 coefficients to describe both components of a 
trajectory, if the database was filling this space, the num- 
ber of points within the sphere would be to oc r 20 (no 
realistic experimental database can "fill" such a volume). 

Conversely, one can define the fractal dimension (or 
Minkowski-Bouligand dimension [13J) p of a set (in partic- 
ular, of a "database" of Laguerre coefficients) by counting 
the number m(r) of pairs of points in the set which have 
a distance smaller than r: 



p(r) 



d log m(r) 
dlogr 



(82) 



One should either smooth m(r) or compute the derivative 
by finite differences over a large enough interval. Note 
that the fractal dimension p is a function of the scale r. 

Imagine that we have a series of data-points (x, y, z), and 
we are investigating whether z can be predicted using x 
and y. Let us imagine that the fractal dimension of the 
set of (x,y) pairs is 2 (the set of (x,y) fills the plane). 
If the fractal dimension of the set of (x,y,z) is equal to 
2, then the set of (x,y,z) is within a surface, and z can 
be predicted using x and y. If the fractal dimension of 
the set of (x,y,z) is equal to 3, the data forms a cloud, 
and x and y are not sufficient to predict z, other hidden 
variables must be at play. These concepts are now applied 
to the study of the database. 




10 10 
Velocity [Re] 



Figure 10: Computing fractal dimensions 



(r = 3 x 10 3 [ilr Re]) one can discern a cloud of dimen- 
sion 4.76. At r = 1.5 x 10 3 [ilr Re] the slope decreases 
to about p = 2, and it is believed that this is the dimen- 
sion of the dataset for a given point along the riser. At 
small scale (r < 1 x 10 3 [ilr Re]), the dimension increases 
again, possibly due to noise in the data, or weaknesses in 
the Laguerre approximation. 



The red curves in Figure 10 are computed by adding the 
sum of squares of the differences between force compo- 
nents (suitably scaled) to the squares of the distances 
between trajectories, and then extracting the square root. 
The four red curves are drawn using the original force 
data, to which Gaussian noise of standard deviation 0, 
10 7 , 10 s and 10 9 [N/m] respectively has been added. 
The standard deviation of the original force is about 
2 x 10 8 [N/m] . The two first red curves are indistinguish- 
able, which seems to indicate that we cannot expect to 
achieve a 10% precision in force predictions. The marked 
difference with curves 3 and 4 shows however that we 
have assets in hand to predict the force. Similar curves 
have been produced with added noise of standard devia- 
tions 1 x 10 7 , 2 x 10 7 ... 10 x 10 7 [N/m], and already at 
2 x 10 7 [N/m] the curve is distinct from the one based 
on the original data. 



7 Dynamic analysis 



7.1 Foreword 



Figure 10 shows the cumulative distribution of the dis- 



tances between trajectories (black curve) computed using 
Equation 79 p is seen to depend on the scale on the 
scale: from afar (r > 2 x 10 [ilr Re]), the slope of 
the curve is zero, hence the dimension is zero: all the 
data are lumped into a point. Zooming into the data set 



Once it is possible to predict hydrodynamic forces on a 
cross section for a given velocity history, the next develop- 
ment is to include the force thus predicted in a dynamic 
time domain simulation. Because the VIV forces intro- 
duce severe non-linearities, a naive connection (where the 
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forces are just added to the right hand side of the sys- 
tem) might lead to slow convergence, or to divergence of 
the Newton-Raphson iterations used at each time step. 
To obtain a proper formulation, it is necessary to jointly 
treat the system of differential equations composed of the 
state equations of the structure, and the differential equa- 
tions ( |42| ) followed by the Laguerre coefficient. However 
in doing so, for each displacement degree of freedom, n 
Laguerre coefficients are added, and it is crucial for ef- 
ficiency to eliminate them before solving a large linear 
system of equations. 

To this effect, in this Section, the following sequence of 
transformations is applied to the differential equations: 

1. The differential equations are first set in incremental 
form (Section |7.3[). 



2. Time discretization by the Newmark-/3 method is in- 



troduced (Section 7.4 1. 
3. The Laguerre coefficients are condensed out of the 



system of equations (Section 7.5 1. 
4. Finite element interpolation is introduced (space dis- 



cretisation), using Gauss quadrature (Section 7.6 1 



This particular sequence leads to a VIV model that is 
implemented at the Gauss point level, and can easily be 
introduced in a general purpose FEM software with stan- 
dard, displacement based beam or cable elements. An- 
other sequence, 1, 4, 2, 3, can be used to obtain either 
a hybrid element, or alternatively, a mixed element which 
would require a specialized solver for optimal efficiency. 
These alternatives are more difficult to integrate into ex- 
isting software working with displacement based elements, 
and are not discussed here. 



7.2 Differential equations 

The dynamic differential equation of a 3D beam subjected 
to VIV loads can be formalized as 



Tdi ( x [bj],i[bj],X[bj],t) 

= ^/Vm-l f d ( T [pb]t) + E di (83) 

where Newton's "dot" notation for a time derivative 
stands for a derivation with respect to unsealed time t, 
as opposed to scaled time t* , and with [17] [5] 

Edi = C L pv (w d i - ±di) 

+ C Q ^pDi\w di - x di \(w di - x di ) 



+ C M -^pD 2 t (w di - x di ) 

7T 2 

+ ^pDiW di 



(84) 



added mass forces, and the Froude-Krylov forces. The 
fourth term introduces the correction discussed in Section 

EH 

If Cx, Cq or Cm are set to values different from zero, 
then it is necessary to substract the correspond values 
from the forces f d i used to train the rotatron. Experience 
shows that the Using Cm = 1, Cq = 1 and Cl = 
contributes to the stability of the dynamic analysis. 



Equation 42 must be scaled to keep only derivatives with 
respect to unsealed time, for the application of Newmark- 



j3 (Section \7A\ 

dribi 
dt* 



so that 



A. 1 fn,i 



= PipTpbi + ni\ m ], (x bi - w bi ) (85) 
MpTpbi + ni\ ms -i (xu - w b i) (86) 



The indices d and b span pairs of directions, orthogonal to 
the cylinder. Indices i and j stand for positions along the 
cylinder, and span a continuous set of values (coordinates 
along the cylinder). Indices I and p refer to the Laguerre 
coefficients of various degrees. Forces ] dl = f d (ripbu) at 
location i only depend on the Laguerre coefficients t p u for 
the same location. At that location, the force component 
in direction d depend on the Laguerre coefficients of all 
degrees b for both directions p. p is the fluid density w d i 
is the acceleration of the undisturbed fluid. jpD^w d i 
stands for the Froude-Krylov forces. Diffraction forces 
are present in the laboratory tests and hence accounted 
for by f d . 



7.3 Incremental form 



The incremental form of Equations 83 and 86 s 



Tdi + k dib jdx b] + c dl bjdx b3 + m di bjdx b j 

= ^ Nm -ifdi + hdipbjdTpbj + E di (87) 

K 1 Cn« + df Ibi) = ( T pbi + dTpbi) 

+ ni\ ms -i (± bi + di bi - w d i) (88) 

(89) 



with 

k d ibj 
Cdibj 



dr di 
dx b j 
dr di 



x bj 



CQpDiSi 



\w 



[/■'] ' 



m-dibj 
hdipbj 



+ (w di - ±di) (wm - x b i) \w[ p ]i - i[p]j 

+CLpv5ij5bd 
dr di it 

-qJT— + CM-^pViOijObd 



- A" 1 



df t 



di 



Nm- 1 



dr. 



pb 



(90) 
(91) 

(92) 



The four terms in the above Morison's equation are the 
linear drag, the quadratic drag, the sum of diffraction and The expression for |^ is presented in Appendix 
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7.4 Time discretization 



7.5 Condensation 



Newmark-/? is a method geared towards 2nd order differ- 
ential equations. Equation|88] however, is only of the first 
order, and this opens two options: we can treat Equation 
as being of the second order in T[ bi , but with the coef- 



ficient of Tim being zero. Alternatively, we can introduce 
the antiderivative Tm of Tm, and treat Equation 88 as 
being of the second order in TJ&j, but with the coefficient 
of Tim being zero. The later option was chosen, based on 
the weak justification that this treats Ti bi and x b j both as 
first derivatives, which seems natural considering Equa- 
tion EH 



Applying Newmark-/? to Equations 87 and 88 in this way 
yields 



Vd, i, 



^dibj 



7 

Pdt 



Cdibj 

7 



1 



Pdt- 



; m dibj 



dx 



Pdt 

-lfdi + Edi — r di 

Cdibjbftj + ITidibjdbj ~ hdipbjbpbj 



- X nI 



(93) 



and 



(3dt 



ni\ ms -idx u + 



1 



X^S, 



pdt 2 " 8 " lp pdt 
ni\ ms -i (x b i - wu) + VipT p u - K 1 
— niX ms -ib bi — [ii p bp bi + 



7 



Hip 



dT, 



Tlbi 



'lb, 
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a bj 


= Wt ib ^ 
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" 2p Xb3 
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= jibj + \ 




T 
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1 

" p~dt TpbJ 


1 . 

+■ Yp T ^ 


°pbj 
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dt 



Xbj 



dt- 



pbj 



For refinement iterations, 



b x b y a T pbj and b T pbj 



pbi 

(94) 

(95) 
(96) 
(97) 
(98) 
are set 



The time discrete equations can be rewritten in a compact 
form: 



3 dibjdXbj 



stdXb 



■'dipbj 
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dT, 



pbj 
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i/ith 
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(103) 
(104) 

(105) 
(106) 



(107) 
(108) 

(109) 

l Tlbi 

(110) 



One can then condense dT^j out of the above system of 
equations: 

-1 1 « 



dTphj — 



dx h j) 



^dibj 



5 'dipbj 



(*%i 



dx 



bj 



dipbj 



-16 
pi 



(111) 



(112) 



Equation |112| is forced into "Newmark" form as 

k 



^dibj 



^dibj 



7 



1 



pdt 



Cdibj 



pdt 2 



m dibj 



dx 



b.) 



v dibj 



= fdi - r di + CdibjKj + m dlbj a bj 



/rfi ~~ ^N m -i fdi + ^rf? 



hdipbj bpbj 



„2 

dipbj 



(* 5 ) 



(113) 
(114) 
(115) 



to zero. Typicaly, 7 =~, /3 = 4. The step dt refers to fc^and both depend on dt, p and 7: The symbol 



aled tir 



dibj' 
^dibj 



Idi 

was chosen to indicate that the matrix is handled 



by the Newmark-/? solver in the same way as a stiffness, 

/\ 1 ■ m 1 o x-v. j x-v. u r however this terms can not be interpreted physically as a 

As usual in the Newmark-p method, the increments tor r r 7 J 

the time derivatives are found from the increment as 





m dxbi - 
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pdt dTpb3 ' 
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dTpbj 


1 rlT 1 

Pdt 2 pb] 


T 

— a pbj 


(102) 



7.6 Spacial discretization 

The consistent discretisation by Galerkin finite elements 
of Equation 113 leads to 



K* = 

nm 



Ndinkrfrtj N b j m 



F n — N din fdi 



(116) 
(117) 
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K* m and F* are typicaly computed by Gauss quadra- 
ture. Note that no space derivative is present in 
so no partial integration or Gauss quadrature with curva- 
ture shape function appears. One can hence simplify the 
expression of the element matrix to 



Test name 


Reynolds number 


Current 


TN2030 


13500 


uniform 


TN2340 


0-16200 


shear 


TN2370 


0-24300 


shear 



,,,,, "din K dib ly bim 



which means "same quadrature as for a mass matrix". 
7.7 Implementation 

In non-linear FEM code, incremental matrices and vectors 
are computed by Gauss quadrature. The Gauss quadra- 
ture involves shape functions, tensors that are local, con- 
tinuous versions of the stiffness, damping and mass ma- 
trices, and the force imbalance vector. For example for 
the drag damping of a beam element, the tensor relates 
a vector which components are increments in velocities in 
three directions, to another vector which components are 
increments in forces per unit length in three directions. 

Within an iteration, the linear solver provides incremen- 
tal nodal positions, velocities and accelerations for the 
model. These are disassembled and provided to the el- 
ements. The elements compute positions, velocities and 
accelerations (and more) in a co-rotated reference system 
at Gauss points. The resulting values are handed to the 
VIV-Gauss point procedure. 



Table 1: Reynolds number in NDP tests 

x io 4 Training set 



-1.5 - 




0.5 1 1.5 
V„ [Re] 



Figure 11: Force vector and Laguerre approximation of 
velocity, for a fraction (1%) of the data used to train the 
rotatron. The blue cross marks the origin (zero velocity 
relative to water) 

8 Results 



The axial velocities are discarded. The procedure scales 

the provided values using Equations [4] and [5] Training 



Having stored the previous approximation of the scaled 
position, the procedure determines the position increment 
dxbj, and then uses Equation | 111 to obtain the Laguerre 
coefficient increment dTp^j. From there, Equations 101 
and 102 are used to compute dr^j and df p bj . The values 
of Tpbj, T p bj and T p bj are updated from previously stored 
values. Tpi,j is then used to evaluate fm and its derivative 
with respect to T p bj. These are scaled back, and Froude- 
Krylov forces are added, leading to k* ub j and j* u . 

The above matrix and vector are padded with zeros to 
indicate zero force in the axial direction and zero torque. 

The condensation of a larger system of time-discretized 
equations introduces some inelegant features compared 
to standard dynamic FEM: the VIV-Gauss point must be 
provided with (3, 7 and dt and a flag showing wether a 
call is made at a step or within a refinement iteration. 

Note that the dynamic FEM computation does not make 



The Norwegian Deepwater Program was a research effort 
in which reduced scale tests were carried out on long, 
flexible riser models, subject to uniform or sheared cur- 
rent [3]- The displacement histories thus aquired at 19 
points along the riser model were later prescribed on short 
stiff cylinders, and the hydrodynamic forces acting on the 
cylinders directly measured [?]. The data from [?] that 
is used in this work consists of the displacements at 19 
points along the NDP riser model, for 3 current profiles 
(Table]!]), so a total of 57 short cylinder runs. For each of 
the 57 runs, 100 instants are randomly selected, yielding 
a training set of the rotatron with 5700 "points". Each 
"point" consistes of two sets of n = 30 Laguerre coeffi- 
cients and the two components of the corresponding force 



(Figure 11 1 



The rotatron was trained using n — 30 Laguerre polyno- 
mials, 200 neurons in the hidden layer, and 50 to 1000 it- 
erations of the conjugate gradient optimization algorithm. 



any use of the recursive filter presented in Section 4.6 this 
filter is used only in the training of the rotatron model. 
In the context of training, the filter had the advantage 
of making refinement iterations unnnecessary. It further 
allowed to avoid using Newmark-/? in a situation with 
prescribed displacement, for which it is not well suited. 



Figure 12 shows how the rotatron predicts the forces for 
the trajectories in the above-mentionned 57 runs of short 
cylinder tests, the comparison of the forces aquired exper- 
imental with the forces predicted using the model. The 
model's ability to predict these forces seems to be good, 
although we lack a good criteria to judge that yet. 
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TN 2030 TN 2370 TIM 2340 



Figure 12: Quality of prediction on the training data set. Velocity (black), training force (red), predicted force (blue). 
All velocities and all forces presented at the same scales 
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x io 4 Time: 875.0 [ji s] Lgr err: 297 [Re] 




-1 - 



5000 10000 15000 20000 

V x [Re] 



Quantity 


Value 




length 


38 


m 


■ 1 ■ . 

outer diameter 


0.027 


m 


EI 


37.2 


Nm 2 


EA 


5.09- 10 5 


N 


mass 


0.933 


kg/m 


tension 


3000 


N 



Table 2: Characteristics of the NDP reduced scale riser 
model 

damping, length and boundary conditions affect the vibra- 
tion and hence the velocity trajectories that appear in the 
vibration. Hence the database used in Section IQl to train 
the rotatron is specialised, not only to a few Reynolds 
numbers for the current velocity, but also to some extend 
to the particular model used in the NDP program. No 
study was carried out in this work on which changes in 
the structure change its way of vibrating to the point were 
the rotatron provides poor for estimations for it. 



Figure 13: Illustration of the approximation process. See 
Section |8.1 on page 15| 



Figure 13 provides a visualisation of the different steps of 



the modelisation process, and is hence a useful diagnostic 
tool. It shows: 



Stippled black line The trajectory for which a force 
prediction is wanted. 

Smooth black line The Laguerre approximation to 
the above trajectory, used to enter the rotatron. 

Stippled black arrow The measured force for the 
above trajectory. 

Smooth black arrow The predicted force for the 
above trajectory. 

Green Neighboring (in the sense of the rotatron distance, 



Equation 79 on page 11 ) trajectories from experi- 



ments, used in the training set (and corresponding 
Laguerre approximation, experimentally measured 
force and predicted force). 

Red Same as the above after rotation and/or mirroring. 



Figure 13 gives an indication of the quality of the Laguerre 
approximation, the adequacy of the training set for the 
trajectory at hand, the presence of contradictions in the 
training set near the trajectory at hand, the quality of 
the fit of the rotatron to the training data and finally the 
quality of the interpolation between training points. 



8.2 Dynamic analysis of a flexible riser 



Hence, in order to test the performance of the present 
VIV model within a dynamic analysis, the simplest case 
was considered: the riser model used in the NDP testing 
program (|3J, characteristics in Table [2]) was modelled. 
The present method was used to analyse the test con- 
dition TN2370 (0 - 24300 Re). Numerical results were 
compared with those obtained experimentaly on the NDP 
model. A laptop, using one core of a dual core processor, 
took 20 to 50 seconds to compute one second of riser 
response. 



In Figures 14 and 15 the horizontal axis is NDP- 



laboratory time, the vertical axis is the unsealed length 
along the riser. The upper subplot shows the response in 
line (IL) with the flow, the lower subplot the cross-flow 
(CF) response. The color codes the displacements, with 
the same color scale used in Figures 14 to 16 



Figures [14] and [15] are for test TN2370 (0 - 24300 Re). 
The dynamic simulations captures the frequency doubling 
between CF and inline, as well as the instationary na- 
ture of the vibrations. Frequency and amplitude are ade- 
quately captured. The 6th mode's dominance of the CF 
vibration is correctly captured. The dynamic analysis as- 
sumes constant tension. By contrast, some small tension 
modulations seem to displace the position of the lower 
vibration node in the test. In test 2370 there is a marked 
tendency for CF vibrations to propagate downwards. In 
the analysis results, CF waves are of a more static char- 
acter. The IL vibration in the analysis occurs at a higher 
mode (11) than in the test (9). This is best seen by 
counting red dots along diagonals, for example starting 
from coordinate 38m and time 17.9s in Figure [14] As a 
consequence, and since, for a given propagation celerity, 
wavelength and period are related, the phase drift be- 
tween IL and CF are of opposite sign in the analysis and 
the test. 



VIV depends not only of current velocities, but on the It proved impossible to reproduce tests 2340 and 2030 
type of slender system they act upon. Tension, stiffness, in the same manner. The analysis quickly ends with IL 
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Figure 14: NDP test results, test 2370 (0 — 24300 Re). The color coding describes the displacement 
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Figure 15: Analysis results, test 2370 (0 - 24300 Re) 
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Figure 16: Analysis results, test 2030 (13500 Re) 



and CF vibrations occurring at the same frequency, and 9.2 Dynamic analysis 



n Figure 16 with in-line motions dominating. 



9 Discussion 



9.1 Force prediction 



Figure 12 shows that given the velocities in the training 
data as input, the model allows to reproduce the forces 
in the training data, based on the velocity in the training 
data. 

The same exercise is carried out with trajectories from 
another test, N2430 (not to be confused with TN2340), 
which was carried out in shear current (ORe to 40500 Re 
). In comparison, the highest current velocity appearing 
in the training set is 24300 Re. In N2430, the forces are 
fairly well predicted at the lower current velocities, while 
the predictions are very poor at higher velocities. This 
illustrates that the present model provides no mechanism 
to "extrapolate" over Reynolds numbers, in contrast to 
other VIV models. One could imagine a future model 
in which some dependency on the Reynolds number is 
encoded, just as symmetries are now encoded in the ro- 
tatron. However, this will not be simple: if for exam- 
ple, viscous forces are to increase quadraticaly with the 
Reynolds number and the inertial forces linearly, how in 
the first place to partition a hydrodynamic force from a 
test in two such components? Would such a scaling work 
at all? 



The comparison of Figure 14 and 15 is encouraging: when 



simulating exactly the system from which trajectories were 
aquired experimentally in the first place, the simulation 
come strickingly close to the experimental results. It is 
understood that a simulation is good if it captures the sta- 
tistical properties of the real response. VIV being chaotic, 
there is no hope to reproduce exactly any given realization 
of the response. 



Figure 16 shows the result of a simulation of a case which 



is also directly represented in the training data. The sim- 
ulation fails, in the sense that IL and CF vibrations are 
simulated to occur at the same frequency. Two explana- 
tions are proposed: 

The quality of the force prediction by the rotatron, illus- 
trated in Figure 12 was declared "satisfactory", but this 
constitutes only an observation that the fitting procedure 
is operating. On what criteria should one judge that the 
fit is adequate? When training was carried out, the av- 
erage norm of the difference between training force and 
predicted force was found to be about 30% of the average 
norm. For a control group composed of data points from 
the same experimental database, but not used in training 
the perceptron, the same ratio was about 40%. These 
number are high, and can be reduced to some extend by 
increasing the number of hidden layers and of training 
iterations, but this was not found to yield better simula- 
tions. Further specialising the perceptron (by training it 
with data from test 2030) did not lead to successful sim- 
ulations. One study which might help to understand the 
observations would be to find the corrective forces that 
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need to be added to the forces predicted by the model, to 
force the simulation to track the motions observed during 
test 2030. This can be achieved using dynamic inverse 
FEM analysis [10] HUE]. 

Even if the above corrective forces were strictly zero, so 
that the model was accurately predicting forces for the 
riser motion from test 2030, this would not be sufficient to 
ensure that the simulation adequately mimics test 2030. 
A given trajectory could still have very different stability 
properties in the physical system and in the simulation. 
It could be that in the physical system, the trajectories 
follows the "bottom of a valley" while model renders it 
as the "crest of a mountain". A measure of stability for 
this is the Lyapunov exponent [9j [20j . Procedures exist 
to compute Lyapunov exponents from experimental data, 
and this should be compared to Lyapunov exponents for 
the simulation. 



9.3 Influence of L, 



7 iR 






One issue that was explored was the adequacy of t w : The 
rotatron was trained with t w = 5-10~ 5 . This corresponds 
roughly to 1/4 of a cross-flow oscillation period in reduced 
scale for test 2370, and to a smaller fraction for other tests 
at lower current velocity. It could be argued that this frac- 
tion becoming too small could be the cause of the failure 



of analyzes at lower velocities (Figure 16 1. The rotatron 



was hence trained again using a suitably increased value 
of t w . This was done twice, once with the full training 
set and increase the number of Laguerre polynomials to 
n = 60, and once with only the experimental data from 
lower currents, and an unchanged number of polynomi- 
als (n = 30). Neither rotatrons allowed to perform a 
successful simulation for the lower current velocities. 



9.4 Tension 





In Figure [14] one can observe modulations of the posi- 
tions of the vibration nodes (in particular on the cross-flow 
graph). One possible explanation for this would be that 
in the physical system, the tension is modulated by the 
vibration. This effect, if present, is not captured by the 
numerical model, which assumes constant tension. The 
method presented here is designed for use in a non-linear 
analysis. However in this research, time was saved by 
using a simpler linear structural model. 



N 2430 



10 Conclusions 



Figure 17: Quality of prediction on a reference data set. A model for the prediction of VIV forces given the his- 
Cf. color code of Figure [12] tory of velocity of a cylindrical cross section relative to 

the undisturbed fluid, has been developed. The model is 
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closely relatied to Wiener-Laguerre filters: the recent his- 
tory of velocity is represented by the coefficients of a La- 
guerre polynomial series. These coeffcients are then used 
to enter a memory-less non-linear interpolation function, 
in this case, a custom made neural network in which some 
relevant symmetry properties were "hard-wired". The neu- 
ral network was trained by using forces and displacements 
obtained in irregular forced motion tests on a short cylin- 
der. 

The proposed model operates in the time domain, making 
it well suited for integration into fully non-linear analyses 
with unsteady currents. It further deals with in-linea and 
cross-flow vibrations as one inseparable issue, which is 
arguably a necessity to improve on existing VIV models. 

The model could provide a "good" reproduction of the 
forces in the training test, as well a "good" prediction of 
forces for "comparable" trajectories. Were the model was 
queried with trajectories very different from those present 
in the training data, the model gave very poor results 
- as can be expected. Due to the limited amount of 
experimental data available, the present model remains 
quite specialized to a limited number of situations. What 
remains unknown at this stage is the size of the training 
set, and of the neural network model, necessary to create 
a model with some pretention of generality. 

In some in dynamic analyses of laboratory tests (NDP 
TN2030 and TN2430) with a long flexible riser, the nu- 
merical solution fell into an unphysical mode of vibration 
with the same frequency for in-line and cross flow vibra- 
tion. On the other hand, for another case (TN2470), 
some unusually fine details were captured by the numeri- 
cal model. 



A Conventions for indexed nota- 
tions 

In the present work, index notations inspired from tensor 
analysis are used. However, the present setting differs 
from tensor analysis in at least three ways: 

First, we assume that we are only operating in Euclidian 
spaces (an not in more general Riemannian manifolds) so 
that orthogonal bases can be used. This makes it unnec- 
essary to distinguish between co- and contravariant bases 
and coordinates. Hence, only lowered indexes appear in 
the present work. Incidentally, it was here assumed that 
the state of the model is a point in a vector space, which is 
not true when finite rotations are present and Riemanian 
geometry should be introduced instead. 

Second, in tensor notations, each index spans the dimen- 
sion of the manifold. In an expression like cry = Cijki £ki 
the indices range from 1 to 3. Following Einstein's con- 
vention, indices k and I are summed over, and the relation 
is valid for any combination of i and j. The fact that the 
equation is valid at each point within a solid is implicit 
in the notation. In the present work we prepare for the 
manipulations of arrays in a computer, involving opera- 
tions that are repeated, for example for various locations 
alomg a riser. If indexes x, y and z were introduced to 
note the position to which the various tensors refer, one 
would tend to write a ijxyz = C l]k i xyz e k i xyz , which vi- 
olates Einstein's convention, because no summation (or 
rather: no integral) is implied over the positions. 

Third, we introduce non-linear functions. These functions 
can combine the values of the coordinates for some in- 
dices, and operate in parallel on the coordinates for other 
indices. 



From this, it is concluded that the concept has merit and Hence the following conventions are used: 
deserves to be pursued, acknowledging that more work is 
needed to arrive to a pratical engineering tool. 

1. By default, where an index appears more than once 
in a combination of products and/or divisions, a sum- 
mation over the index is implied. If that index has 
a continuous range, then the "sum" is an integration 
over the range. Point 2, 3 and 4 specify exceptions 
to this rule. 
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2. Point 1 notwithstanding, if an equation is preceded 
by the symbol V, followed by a list of indexes, then 
the listed indexes are not summed over. 

3. Point 1 notwithstanding, if within an equation, there 
is a combination of products and/or divisions within 
which an index appears only once, then no sum over 
that index is carried out in the whole equation. (In 
any other situation than a simple term in the left 
hand side, readability should be improved by using 
the symbol V . ) 

4. Point 1 notwithstanding, if an index appears within 
an input to a function, and the output of the function 



is multiplied or divided by one or several terms that 
have the same index, then no sum within the input 
to the function is carried out on that index. 

If an index of an argument to a function is within 
brackets, then the whole range of index values is used 
as input to one function evaluation. For example, 
°i (y[j]k) refers to the evaluation a multiple locations 
(k) of a vector-valued (i) function of a vector (j). 

When the output of the function is shorthanded with- 
out explicitly writing its input, then the indices of the 
input that are not within bracket are added to the 
indices of the function. For example Cj (yu]k) can 
be shorthanded <T ik . 

Derivatives of a function are noted with only the 
bracketed indices of the input appearing under the 
fraction: S^* . To refer to the value of that deriva- 

tive for input k, one writes . 



B Rotatron gradients 



with 

d<T ik r 



dU, 



-5, 



kr 



e l Uikn \y[j]kn 



\OLk-l 



log \V[j]kr, 



{ I I a fc -1 \ * 

{\y[j]kn\ +i) 



(126) 



C Inverse of s5 



The inverse of s5 (Equation 109 1, where s5 is of the form 

s5 y = aT i:j + pSij (127) 

with 

f/ " = 1 ; ; (128) 



= 



3 < i 
3 > i 



can be verified to be lower triangular banded, with terms 
on diagonal i equal to 
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e{2...n} 



(129) 
(130) 



The derivative of the force predicted by the rotatron, with 
respect to the Laguerre coefficients is needed in Section 
[7] With references to Equations 64 to 69 that describe 
the rotatron, we can write 



dfin 



= Vi 



dfin d&ikn dyjkn 

da k dyj di\ 

9(Jikr 



ki 



(119) 
(120) 



with 

d(Jikn 

dy 3 



I I 3 !\ l Q k 
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1)' 



I a k 



&k Uikn Ujkn \U[j]k7i \ 

+(-l) s 'iy^ kn y^kn {\y\j]kn\ a " + l ) } ( 121 ) 

Here index i ranges over two values (for two directions 
orthogonal to the cylinder), and -^i is the other direction 
than i. 

The gradients of the rotatron with respect to its coeffi- 
cients are also needed in order to compute the gradient 
of the target function with respect to the parameters V k , 

U k and M M . 
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